# 
# data <- readRDS('data_sim_p_1000.rds')
# for(x in names(data)) assign(x, data[[x]], envir = globalenv())
# U_full <- data$U
# rm(data)
# 
# S.data.time <- S.data$obs
# S.data.time.log.sum <- sum(log(S.data.time))
# Exemple d'optimalité entre les 3 fct MH, MH_high_dim et MH_high_dim_para sur le modèle non lineaire mixte
m <- function(t, phi1, phi2, phi3) (phi1  )/(1+exp((phi2-t)/phi3))
#=======================================#
p <- 1000
parameter <- list(sigma2 = .05^2,
                  mu = c(0.9,90,5),
                  omega2 = c(0.005, 40, 1),
                  bara = 90,
                  barb = 30,
                  baralpha = 0.5,
                  beta = rep(0,p))
parameter$beta[1:4] <- c(-.8, -.2 , .3 , .9)
#=======================================#
t <- seq(60,120, length.out = 10) #time values
set.seed(123)

G <- 40 ; ng = 4
link = m

dt <- create_JLS_HD_data(G, ng, t, m, link, parameter)

var.true <- dt$var.true
a <- var.true$a ; var.true$a <- NULL #a fixé (et retiré des variables latentes)
S.data <- dt$survival
U_full <- dt$U

Y <- do.call(get_obs, var.true) + rnorm(n, 0, sqrt(parameter$sigma2))
S.data.time <- S.data$obs
S.data.time.log.sum <- sum(log(S.data.time))
model <- SAEM_model( 
  function(sigma2, ...) -n/(2*sigma2),
  function(phi1, phi2, phi3, ...) mean((Y - get_obs(phi1, phi2, phi3) )^2 ), 'sigma2',
  
  # === Variable Latente === #
  latent_vars = list(
    # === Non linear model === #
    latent_variable('phi', dim = G, size = 3, prior = list(mean = 'mu', variance = 'omega2'),
                    add_on = c('zeta(phi1 = phi1, phi2 = phi2, phi3 = phi3, ...)' )),
    
    # === S.data model === #
    latent_variable('b', prior = list(mean = 'barb', variance.hyper = 'sigma2_b'),
                    add_on = c('zeta(b = b, ...) +',
                               'sum(h$eval(b = b, ..., i = c(1,2)))' )),
    latent_variable('alpha', prior = list(mean = 'baralpha', variance.hyper = 'sigma2_alpha'),
                    add_on = c('zeta(alpha = alpha, ...) +',
                               'alpha*h$eval(alpha = alpha,..., i = 3)'))
  ),

  # === Paramètre de regression === #
  regression.parameter = list(
    regression_parameter('beta', 1, function(...) SPGD(5, theta0 = beta,
                                                      step = 0.05, lambda = 1/sqrt(N),
                                                      normalized.grad = T,
                                                      zeta.der.B, N, zeta.B, 
                                                      Z$alpha,  Z$phi1, Z$phi2, Z$phi3,Z$b) )
  )
)
# ---  Initialisation des paramètres --- #
parameter0 <- parameter %>% sapply(function(x) x* runif(1, 1.1,1.4))
#===============================================#
load.SAEM(model)
U <-  U_full
oracle <- maximisation(1, do.call(S$eval, var.true), parameter, var.true)
#==============================================================================#

init.options <- list(x0 = list(phi = c(1,80,4), b = parameter0$barb, alpha = parameter0$baralpha), 
                     sd = list(phi = c(.05, 1.5, .5), b = 1, alpha = .1) )

SAEM.options <- list(niter = 300, sim.iter = 5, burnin = 190, 
                adptative.sd = 0.6)

saem

p <- 4
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 48min 37sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 0.5000
Estimated value 0.0025 0.9033 89.9108 5.0104 0.0040 35.2265 0.5900 29.6106 0.7293
Rrmse 0.0022 0.0001 0.0005 0.0005 0.0105 0.0192 0.1507 0.0130 0.4585 
p <- 20
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 48min 51sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 0.5000
Estimated value 0.0025 0.9033 89.9168 5.0148 0.0040 35.1852 0.5881 28.5855 0.8915
Rrmse 0.0025 0.0002 0.0005 0.0014 0.0207 0.0204 0.1535 0.0471 0.7830 
p <- 100
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 54min 04sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 0.5000
Estimated value 0.0025 0.9035 89.9555 5.0277 0.0039 35.2423 0.6457 44.0060 2.1571
Rrmse 0.0006 0.0003 0.0000 0.0040 0.0095 0.0188 0.0706 0.4669 3.3142 
p <- 500
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 55min 16sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.000 0.5000
Estimated value 0.0025 0.9034 89.9250 5.0230 0.0039 35.5875 0.5844 35.969 1.3529
Rrmse 0.0017 0.0003 0.0004 0.0030 0.0007 0.0092 0.1588 0.199 1.7059 
p <- 1000
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 54min 29sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 0.5000
Estimated value 0.0025 0.9039 89.9224 5.0250 0.0039 35.3015 0.5846 30.6821 0.6981
Rrmse 0.0018 0.0008 0.0004 0.0034 0.0026 0.0172 0.1585 0.0227 0.3963